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Abstract. We propose a fast multi-orbital impurity solver for the dynamical mean 
field theory (DMFT). Our DMFT solver is based on the equations of motion (EOM) 
for local Green's fimctions and constructed by generalizing from the single-orbital case 
to the multi-orbital case with inclusion of the inter-orbital hybridizations and applying 
a mean field approximation to the inter-orbital Coulomb interactions. The two-orbital 
Hubbard model is studied using this impurity solver within a large range of parameters. 
The Mott metal-insulator transition and the quasiparticle peak arc well described. A 
comparison of the EOM method with the QMC method is made for the two-orbital 
Hubbard model and a good agreement is obtained. The developed method hence holds 
promise as a fast DMFT impurity solver in studies of strongly correlated systems. 



PACS numbers: 71.27.-|-a, 71.30.-Fh, 71.10.Fd, 71.10.-w, 71.15.-m 



Fast multi-orbital equation of motion impurity solver for DMFT 



2 



1. Introduction 

Since the introduction of the dynamical mean field theory (DMFT) [H |2] in the past 
decade, it has been continuously developed and improved and has become a powerful 
tool to investigate strongly correlated systems involving d and / electrons. DMFT 
greatly helps to understand, e.g., the Mott metal-insulator transition as well as other 
exotic properties in strongly correlated electron systems. The basic idea of the DMFT 
is to map interactions between different sites to interactions between an impurity and a 
self-consistently determined bath |3]. Then, a lattice problem is transformed to solving 
a single impurity problem along with a DMFT self-consistency condition. Thus, the 
central part of DMFT calculation is to solve the impurity self-energy efficiently. This 
is a numerical problem, to which a considerable amount of effort has been devoted to 
advance and improve it during the last decade (see, e.g., [1] for a review). 

Frequently used impurity solvers are the quantum Monte Carlo (QMC) method 
[5], the exact diagonalization (ED) method [HI E], numerical renormalization group 
(NRG) method [H [9], density matrix renormalization group (DMRG) method [10], 
and recent continuous time quantum Monte Carlo (CTQMC) method [lU |T2]. But 
all of these methods are either computationally expensive or work only in a limited 
parameter region. Therefore, novel impurity solvers that are both fast and reliable 
are demanded, particularly when calculations for real materials are performed, merging 
density functional theory (DFT) plus DMFT approaches. The equations of motion 
(EOM) method and associated methods have been studied a lot in the past half century 
[la [H Iia [la [HI [111 [la EQl Ell [221 been considered to be a 

good candidate for a fast impurity solver. In [27] and [2H], an infinite U and a finite U 
impurity solver were successfully developed for the single-orbital case based on equations 
of motion of the Green's functions. Another recent work by Zhuang et al. exploring 
multi-orbital impurity solvers is based on the Gutzwiller variational approach |29j . 

The EOM impurity solver has previously been developed for the single orbital 
case [2Sl ED]- However, in calculations of real materials involving strongly correlated 
electrons, one usually has orbital degrees of freedom, i.e., near the Fermi surface there 
exists more than one orbital. For example, when studying transition metal oxides and 
lanthanides, the valence d (and /) electrons usually occupy several, up to 5 (and 7) 
orbitals. These orbitals have similar energies and, consequently, all of these orbitals 
have to be considered in the manifold of localized states. Therefore, it is essential to 
develop the impurity solver from the single-orbital case to a multi-orbital EOM (MO- 
EOM) impurity solver that is applicable in concrete DFT+DMFT calculations for real 
materials. 

For QMC, ED and CTQMC methods multi-orbital solvers have been developed. 
However, according to the earlier investigation for the single orbital case [28l 130] . the 
EOM impurity solver is computationally much faster compared to QMC and NRG 
methods, because the EOM method does not need matrix calculations, but solves the 
integral equations. As in general calculations of multi-orbital systems will be much 
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more time consuming than those of single-orbital systems, due to the increase of the 
dimension of matrix, the MO-EOM impurity solver will hopefully save more cpu time. 
Moreover, it has the advantage that it works directly on real frequencies and can hence 
be employed for an arbitrary finite temperature. 

Our aim is to obtain a fast, reliable, and general multi-orbital impurity solver with 
EOM method. Comparing to the single-orbital system with SU(^A^ symmetry which 
can approximately deal with the A^-orbital systems with degenerated energies, in the 
multi-orbital impurity solver all the interactions of the localized electrons should take 
into account the orbital degrees of freedom and all the orbitals may also variate in 
their physical parameters, e.g., band widths, intra-orbital and inter-orbital Coulomb 
interaction strengths, and so on. With more flexibility in the interactions, a wider range 
of physically relevant cases and more complex systems can be treated, and the method 
can provide more interesting information. 

Furthermore, in all previous studies of the multi-orbital Hubbard model, for the 
hopping term between different sites, the electrons are all treated such that they 
only hop to the orbitals with identical orbital indices on other sites, thus the inter- 
site inter-orbital hoppings are neglected. We think that these inter-site inter-orbital 
hoppings are important because only when all inter-site hoppings combined in one are 
included, can the model accurately reflect a realistic physical image of multi-orbital 
systems. Especially when the orbitals have different band widths the neglect of inter- 
site inter-orbital hopping introduces a man-made difference. Therefore, we have tried 
to incorporate this inter-site inter-orbital hopping effect in our treatment to study the 
orbital selective Mott transition (OSMT) in a fully self-consistent way. 

The paper is organized as follows: in section [2] we describe our MO-EOM method 
and introduce the equations, decoupling, and approximations. In section [3] we present 
model calculation results for the MO-EOM by applying it to two-orbital Hubbard model 
and we compare the computed results with those obtained with the QMC method. Our 
findings are summarized in section |4j 



2. Description of MO-EOM method 

We start from the Hamiltonian of a modified multi-orbital Hubbard model. It is given 

as 

'H = - tijlm filer fjma + Uuhn^fhui + Ulm 

ijlma,iy^j il ilmcrcr' ,l<m 

~^ ^ ^ i.^lma f ima f ila ~^ ^Ima f ila f ima) ^ ('^) 
ilma,l<m 

where 2, j are site indices, /, m are orbital indices, and a, a' are spin indices, respectively. 
The first summation, which sums over two orbital indices and which differs from the 
usually studied form of the multi-orbital Hubbard model [31], is the hopping of electrons 
between different sites. The second (third) summation term is the on-site intra-orbital 
(inter-orbital) Coulomb interaction term, where the intra-orbital Coulomb interaction 
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strength Uu and the inter-orbital Coulomb strength Uimaa' are free parameters that can 
be set all identical or all different in order to simulate different kinds of physical systems. 
The last one is the on-site hopping of electrons between different orbitals which should 
imperatively be taken into account when multi-orbital effects are studied. Here V/^^ 
and V/j^^ are the on-site inter-orbital hopping parameters. 

In the dynamical mean field theory, the many body interactions between different 
sites are mapped to interactions between an impurity and a bath. The most frequently 
used impurity model is the single impurity Anderson model (SIAM). The Hamiltonian 
of the multi-orbital SIAM is 

kla la I 

~l~ ^ ^ Ulmaa'T^laf^ina' 
Imacr' ,l<m 

+ Yl (^lla4kafl- + Vlkaflcika) 
Ika 

+ 5^ (yimafmaflf^ + ^Imaflafma)- (2) 

lma,l<m 

Here hi^j = fj^fi^ is the occupation number for localized electrons with spin a in the 
l-th orbital. The first summation term is the energy of conduction electrons, where the 
electrons in different orbitals are labeled with orbital index I. The second summation 
term is the energy of localized electrons and Sfia- is the orbital level for spin a in the 
l-th orbital. The third summation term is the on-site intra-orbital Coulomb interaction 
term. The fourth summation term is the on-site inter-orbital Coulomb interactions 
between electrons of the l-th orbital and m-th orbital. The fifth summation term is the 
hybridization between the localized electrons and the baths. The sixth summation has 
the same meaning as in ([T]). 

The temperature dependent retarded two-time Green's function in the Zubarev 
notation ^32J is given by 

GAB{t,t') = <^A{t)-B{t') > 

= -iQ{t-t'){[A{t),B{t')]^), (3) 

where A{t) and B{t') are Heisenberg operators, and Q{t — t') is the Heavyside function. 
Here we use the Fourier transform of the Green's function in u space, it satisfies the 
equations of motion 

a;< A;fi>= ([A, 5]+)+ < [A, H^^p]; 5 >, (4) 

w < A; 5 >= {[A, B]+)+ < A- [Himp, B] > . (5) 

In the following calculations of the equations of motion Q is applied. 

For a multi-orbital system, if we define the anti-commutation relation for the 
operators as follows, 

[fla, fma']+ = 0, [//^, /^^/]+ = 0, 

[//o-5 fma']+ = ^Im^aa'i [fla, /mcr']+ ~ ^Im^aa'i 
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we can obtain the first two equations of motion for single particle Green's functions, 

{UJ + H- Sfma) < fnia] fLa > = 1 + Umm < flma'fma] flicj > 

+ ^ {Ulmaa < fli^fma] fLa > 

+ Ulma'a < fli^' fma] fLa > ) 

it 

- 5^ C</;<x;/L>), (6) 

(u + — ekma) ^ Cmka', fma ^ = "^'mmfccr ^ /mo-! /mcr ^ + '^mlka ^ /io"! /mcr ^ •(''') 

Here /i is the chemical potential, a and a' are spin indices where by notation we 
imply a' 7^ a. This notation is employed in the following derivation, too. On the right 
hand side (RHS) of (|6|, the physical meaning of the second term reflects the fluctuation 
of spin a in m-th orbital accompanied by spin a' in m-th orbital, or we can call it the 
fluctuation of spin a when spin a' exists. The third term shows the fluctuation of spin a 
in m-th orbital when spin a in /-th orbital exists, and similarly the fourth term is for the 
case when cr' in /-th orbital exists. The first term (i.e., 1) reflects the existence of spin a 
in the m-th orbital itself. To describe the multi-orbital effect, we have to consider each 
term for the case that other spin-orbital channels exist at the same time. Moreover, ([T]) 
describes the hopping of electrons from the bath to the impurity here. The electrons 
can hop back to any orbital of the impurity. Therefore, from Eqs. ([6]) and ([T]), we can 
imagine that the electrons in the bath which came from the m-th orbital of the impurity 
hop to all orbitals of the impurity, where the V^^f,^ and V^^f^^ are the hybridizations of 
electrons in the bath which came from m-th orbital to m-th and Z-th orbital accordingly. 
There exists the relation that V^,^^ = V^^^^ + J2i,i^mKmka- This constitutes indirect 
hybridizations between different orbitals through the bath, which is generated from the 
bath-impurity hybridization terms of ^ and hence, it reflects the inter-site hopping 
matrix tiji^ in the Hamiltonian of the Hubbard model, ([T]). This is the reason that we 
introduce the inter-site hopping matrix for multi-orbitals in the Hubbard model. But 
one can recognize that V^^^^. and Vf^;,^ will not appear in the impurity Hamiltonian 
(|2| because for indirect inter-orbital hybridizations each orbital of the impurity can only 
"see" the bath. These diagonal and off-diagonal hybridization parameters can only be 
obtained from the mapping of the inter-site hopping on the lattice model. 

In the EOM method when the equations of motion are derived, a decoupling scheme 
is implemented to make the equations closed and solvable. However, this decoupling 
scheme is not only a simplifying numerical technique. We emphasize here that it is 
based on different grades of realistic physical assumptions. In the EOM method the 
interactions included are not simply considered as the interactions shown in the model 
Hamiltonian, but they also depend on the included equations of motion. In calculating 
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the equations of motion, more and more higher order Green's functions will appear in 
higher and even higher order equations of motion. These higher Green's functions are 
associated with some kind of higher order interactions. From the order of the equations 
of motion in which a Green's function first appears, one can define the order of the 
Green's functions and approximately evaluate the contribution of the related interactions 
because this 'order' approximately describes the weight of an interaction. For two 
Green's functions with different order, the lower order Green's function associated with 
an interaction will have a larger weight and hence give more contribution than the 
higher order Green's function associated with the interaction. This is approximately 
valid because the contribution will also depend on the interaction strengths (e.g., the 
values of U, J, V, and so on), and not only on the order. Once one higher order 
Green's function appears, one new higher equation of motion can be derived, while for 
one new equation of motion more higher order Green's functions will appear. This is an 
infinite procedure. The decoupling scheme implies to make a decision where to truncate 
this procedure and only the interactions before this truncation will be exactly taken 
into account and all the interactions higher than this truncation will be approximated 
with lower order Green's functions and relative correlation functions. Therefore, the 
decoupling scheme decides how many interactions are fully included. Understandably, 
a higher order decoupling scheme will take more interactions exactly into account and 
will be more accurate. This is the physical meaning of the decoupling scheme, and it 
is also the reason that we call one interaction a higher order interaction than another 
interaction, by just comparing the orders of the two Green's functions associated with 
the interaction. Thus inclusion of higher-order interactions is a controllable procedure. 
Therefore, the EOM method is a method that can infinitely approach the exact physics, 
whose accuracy depends on the order of the equations of motion or the number of 
interactions that have been exactly included. 

To illustrate this we adopt the single orbital system [28] as an example to show 
the physical implications of the decoupling scheme. For the lowest order decoupling 
which corresponds to the mean field approximation for the Coulomb interactions, all the 
electrons in one orbital are considered moving in the mean field of other electrons, where 
only the on-site Coulomb interactions and the single-electron hybridization interaction 
f^Cka are fully taken into account. Consequently ^ Uafa] fa ^ can be decoupled as 
^o- ^ /o-; /a ^ and only two equations of motion need to be derived, 

+ > = 1 + f/ < ha' fa, fl > +Vka < Cka] fl >, 

(W + /i - < fa] fl^ = < fo] fl >, 

where <^ fia'fa', fl ^ and ^ Cka] fl ^ are the same order Green's functions because 
they appear simultaneously, and are corresponding to the two most basic interactions 
mentioned above. Therefore, the single particle Green's function should be 

« /] »= — ^ . (8) 

CJ + yU — £/ — A — t/ 

which gives only an energy shift of the band along with the Coulomb interaction. 
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For the equivalent of the Hubbard-I approximation which is the second order EOM 
decouphng scheme (here we call this the equivalent Hubbard-I approximation because 
the Hubbard-I approximation is originally for the Hubbard model but not for SIAM), 
one more four-operator single-electron hybridization interaction ficr'flcka is exactly taken 
into account, which relates to the Green's function ^ n^'Cka] fl This interaction can 
be viewed as one next higher order interaction than the Coulomb interactions because 
the Coulomb interaction first appear in the first order equations of motion and this 
interaction first appears in the second order equations of motion. In this case four 
equations of motion will be needed, and finally lead to 

« fl »= . + — . (9) 

cu + fi — Sf — A u + fi — ef — A — U 
which gives the lower and upper Hubbard bands in the atomic limit. However, this 
approximation can not explain the Kondo physics because it can not produce the Kondo 
peak at the Fermi level. Consequently, decoupling schemes beyond the equivalent of the 
Hubbard-I approximation are needed in order to describe the metallic states and explain 
the Mott metal-insulator transition. 

In actual materials calculations, the choice of the decoupling scheme is mainly 
decided by considering two factors: the accuracy and the cpu time consumption. The 
accuracy decides whether the exactly included interactions can be sufficient to describe a 
system. For example, for a complete insulator, the equivalent Hubbard-I approximation 
is already enough and will save a lot of cpu time over the use of the higher order 
decouphng schemes. But for metal and Kondo physics, the decoupling scheme beyond 
the Hubbard-I has to be employed. Moreover, by comparing different order of decoupling 
schemes, it is easy to examine the contribution of the involved interactions, which may 
provide understanding of the underlying many-body physics. 

Now we return to our studied multi-orbital problem. For simplicity, here we first 
treat the on-site inter-orbital Coulomb interaction with a mean field approximation 
which is valid when the inter-orbital fiuctuations are weak in the system, i.e., we use 
the corresponding decoupling 

f^lafmai fma ^ ~ ^'o" ^ fmai fma 
< nia'fmcr; fLa > ~ ^l<r' < fmcT] flia > • (10) 

With this decoupling the on-site inter-orbital fiuctuations are neglected, which may be 
a loss of some interesting information. These on-site inter-orbital fiuctuations will be 
taken into account in a forthcoming work. 

According to ([6]) and ([T]), now the hybridization function A in single orbital case 
will change to the hybridization function Amo- for the m-th orbital. 



T/ 1 )* 

f mfccr ^mmka 



^ U -\- fi — Skma 
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i^lma =2. \ ~ ~, 7t ' ^^^) 

where the ^mma {^ima) are the on-site indirect identical orbital (inter-orbital) 
hybridizations, which relate to the inter-site diagonal (off-diagonal) hopping terms. 
First we shall not discuss the effect of Amma and Aima, and use only the symbol of 
total hybridization Ama in solving the equations of motion. We will then discuss the 
parts of Ama later in this section along with the DMFT self-consistency conditions. 
Thus ([7]) will be replaced with, 

J2 ^^ka < Cmka] fL >= ^rr.a < fma] fL > ■ (12) 
k 

Moreover, because we will mainly concentrate on the on-site Coulomb interactions 
and inter-site hopping effects in ([T]), we first neglect the on-site direct inter-orbital 
hopping terms, which will be studied in a forthcoming work. 

Then we can obtain 

(a; + /X — Efma — ^ {Ulmaa^la + Ulmaa'^la') — Ama) X 
l,l^m 

^ fma'i fma ^= ^ + ^ Uma' fma'i fma ^5 (13) 



where 



ni^ =< ni^ >= j dufiuj - fx)lm<^ fi^-jl:^ (14) 



is occupation number in the l-th orbital and f{uj) is the Fermi distribution function. 

We can see that the existence of electrons in other orbitals has promoted the initial 
orbital level in single orbital case Efma to an orbital level in multi-orbital case, 

Efma = £fma + {Ulmaa^la + Ulmaa'^la')- (15) 

/, Ij^m 

In real calculations, due to the partial filling and the fact that this inter-orbital Coulomb 
interaction should act on the "charge center" of the m-th orbital, this will be modified 
as 

Efma = £fma + (1 " ^ma') X 

{Ulmaaniaflma + Ul 

maa'^la'^ma ), (16) 

I, l^m 

an initial value of which can be extracted directly from density functional theory 
calculations, in contrast to Efma, when we study a multi-orbital system. Afterwards, 
the orbital levels will be corrected iteratively by the occupations obtained in the DMFT 
calculation in each iteration. Moreover, we note that, for a system with fixed filling, the 
Fermi level should also shift accordingly. 

Due to the existence of charges in other orbitals, the energy to add one electron 
into the m-th orbital has also changed from Umm to an effective value 

f/eT = Umm + {Ul nia + Ulmaa'^la')- (17) 

l,l^m 



Fast multi-orbital equation of motion impurity solver for DMFT 9 
For higher order Green's functions, we obtain the equations of motion 



{UJ + ^ — efma)Gnf — ''^ma' + Umm ^ f^ma' fmai fma ^ 

fma', fmc ^ -^Uima'a ^ nia'TT-ma' fma', fmc ^ ) 

I 

+ ( ~ ^mka' ^ (^Inka' fma' fma] fma ^ '^^mka ^ f^ma' Cmka'i fma ^ 

k 

+ Vmka' < /L' Cmfca' fma]fL^)^ (18) 

+ /X — = ^mmka ^ f^ma' fma] fma ^ + O^mk'a' ^ fmcr'^mk'a'Cmka] fma ^ 

fc' 

~ Vmk'a' ^ Cmk'a'fma'Cmka] fma ^ ) "I" ^Imka ^ f^ma' fla] fma 

(19) 

(a; + + - £fema' - ^fma)G'%f = {fLa'^rnka') + V^mifca' ^ fi-ma' fma] fma > 

+ ^ (K^feV ^ fma'^rnka'Cmk'a] fma ^ 

fc' 

~ Vmk'a' ^ C^mk'a''^rak(T' fma] fma ^ ) 

+ I] VLfca' < fla'fl'^'fma] fL >, (20) 

(a; + + - £/mcT' - ^fma)G'^ff = {Cmka'fma') + Umm < Cmka' fma' fma] fma ^ 

"I" 2 ^ ] {Ulmaa f^la^^mka' fma' fma] fma ^ 

+ Ulma'a ^ f^la'C^^j^^i fma' fma] fma ) 

VjTimfco-' ^ f^ma'fmat fma 

+ (Knfe'o-' ^ Cmka' ^mk' a' fma] fma ^ 
fc' 

+ ^mk'a ^ (^mka'fi^a'Cmk'a] fma ) 

- ^ V/mfc(T' < fja'fma'fma] fma (21) 

where we have used the following abbreviations on the left hand side (LHS) of the 
equations: 

^n/ ~ 'f^ma' fma] fma 

^nc ~ 'f^ma'C-mka] fma 

^/c/ ~ ^ fiia' Cmka' fma] fma ^) 

^c// ~ (^mka' fma' fma] fma ^ ■ 

For the three-particle Green's functions, we introduce according to the mean field 
approximation for inter-band Coulomb interactions, 

< rilahma'fma] fLa > ~ ^Z<r < flma'fma] fL (22) 
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< nla'hma'fr,^a^ fl^ > ~ Ul^, < flma'fma^ fla 



(23) 



the effect of this is to increase the impurity position in the same way as in (15) when it 
moves to the LHS in Eqs. (18)-(21) as pre-existed charge in other orbitals. For the same 
reason the Coulomb interaction strength Umm will change to Moreover, if there is 

no external field existing, the energy levels for different spins in the same orbital should 
have Efma' = Ej^a- 

For the Green's functions involving an odd number of operators of the m-th orbital, 
we follow our previous assumption in (11) and treat them as zero. For the Green's 
functions involving only operators of the m-th orbital, we use the decouphng applied in 
[28]. 

Solving now the closed set of equations, we obtain the single particle Green's 
functions 



Tfmm 

^ ^ IP. rrnim A A . A I "'ma' T -( 1 



„ „-j. i^+f^—EfmcT — U^Jg^—AmcT—^ma-l—^rniT \ / fOA\ 



where 



U + fl- Efraa - ^ma - , , , „ p Jjm,r.'" A A A (^"''^ ' h + h 



V* V 



^ W + /i + tjfma' — ^fma — ^kma' 



ma / J 



mka' mmka' 



— U + fl + Skma' — Efma' — Efma — ^If"^ 

I ^mka'^lmka' fl^''i fma' 

kkm ^ + ^ + - - - «fma';fL>> 

_ ^ ^mka' { fma' ^mka' ) ^mka' {(^mka' fma' ) ^ ^27^ 

OJ + 11 + Efma' — Efma — ^kma CO + H + ^kma — Efma' — Efma — U^""^ 

_ — '' ^mka''^mka'{'^mk'a''^mka') , ^mka'^mka' i^mka' ^mk' a') /t)Q\ 



kk 



^ + /i + Efma' — Efma — ^kma W + /i + Skma — Efma' — Efma — Ul 



— Umm + 2 {Ulmaania + Ulmaa'^la')- 

l,l^m 

and furthermore we will use the Hermitian relations 

ifma'^mka') ~ i^mka' fma') ^ i^nik' a'^mka') ~ i^mka'^mk'a') ■ 



In (24), if h, I2 are taken as zero, it turns to a multi-orbital Hubbard-I method 
with the mean field approximation to the inter-orbital Coulomb interactions. , which 
will be reliable at an atomic limit, i.e., for completely insulating states. 
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Let us discuss now on the hybridizations. If we neglect the hybridization between 
different orbitals (i.e., neglect the inter-site inter-orbital hopping in ([T]), and then the 
model Hamiltonian of ([T| will return to the Hamiltonian in [31]), the above formula 
can be simplified by setting Vf^^^ = 0. In this sense, the approximation is similar to 
the single orbital case with arbitrary spin-orbital degeneracy N, where the degenerated 
electrons only hybridize with electrons having identical spin and orbital indices in the 
bath. The main difference between the A^-fold single orbital case and the above multi- 
orbital impurity solver is that the latter one is generalized to A^/2-orbital system with 
different orbital levels, band widths and varied inter-orbital and intra-orbital Coulomb 
interaction strengths for different orbitals. For this case, the DMFT self-consistency 
condition on the Bethe lattice should be 

^ma ~ ^ fmcri fma (29) 

just like the DMFT self-consistency condition for a single orbital Hubbard model with 
SU(A^) symmetry. 

Now let us study the equations in a further step as we mentioned above. If we take 
into account the indirect hybridizations between different orbitals, we will obtain two 
new equations of motion as follows, 

{u + fi- Efi^) < fi,; fl, > = flJ+) + Un < hi,,fi^- /,L > 

+ Uiv„„> < ni,^,fi„] fl^ > ) 
+ 5^ ^fca < q,,; fl^ >, (30) 

k 

{U + fj.- £kla) < C^fca! fLa > = '^Imka < fma'^ flia > +^Uka < //a! flna '^l^m 

+ Y. ^ fl'-^ f^rna > I'^l ■ (31) 

I, Vi^m 

Within the mean field approximation, we can then obtain after charge correction (see 



^ /icr! fma ^ 

/L] + )+ ( Efe Xl^-Zl ) ^ •^«'<^' •^"'^ ^ +(^+M-i?„,-C/"„-A,,-A,,i-A,<,^y , 

= u 1 (^) 

UJ + ^l- hfi, - Ai, - ^^^^E^^^_un^_^^^_^^^^_K,^^la l^k ZJTJI^. 

where we have taken the two correlation functions as 

ifl'Clka') = 0, (33) 

{4k'a'Cika') = - ^ / duj'fiuj' - /i) Im . (34) 

It is a general way to exactly solve the hybridization functions A^a by including 
(32) which will make the equations fully closed for the impurity solver. However the set 



of equations (32) only relate to the hybridization functions. If one has another method 
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to obtain the hybridization functions, the set of equations from (24) to (28) will be closed 
where A^ai and Ama can be obtained from A^a with interpolation, and the correlation 
functions in (27) and (28) are calculated with the obtained hybridization functions 
and single particle Green's functions through the spectral theorem. Fortunately, we 
are now studying a many body system using DMFT. Therefore, the hybridization 
functions Amu do not need to be exactly solved in the impurity solver because in the 
hybridization functions in DMFT will be eventually determined by the DMFT self- 
consistency conditions in each DMFT iteration. The exact solution of the hybridization 
functions will however be needed when only the single impurity Anderson model is 
studied. For DMFT, even if we obtained the exact hybridization functions of the 
impurity, these are only used in the first iteration (guess of the initial GFs) and 
then replaced by new hybridization functions determined by DMFT self-consistency 
conditions from the second DMFT iteration, while using an exact hybridization function 
or an approximate hybridization function in (24) in the first DMFT iteration makes 
little difference for the final result. Moreover, for orthogonal orbitals, the off-diagonal 
single particle GFs need not be calculated. Principally in DMFT for a general lattice, 
the hybridization function should be determined self-consistently together with the 
Green's function and self-energy by summing over k throughout the whole Brillouin zone. 
However, for the Bethe lattice, it will be much easier because it is simply determined by 
the inter-site hopping. Hence, we use the Bethe lattice here as an example. Assuming 
the probability of electrons outgoing from one orbital to a neighbour site is identical to 
the probability of electrons coming from the neighbour site into the orbital with identical 
orbital index. This is reasonable for an equilibrium system. Then, for the Bethe lattice, 
the DMFT self-consistency condition approximately is 

^2 ^2 
Ama = Amma + ^imcr = tm72~ ^ fma'j fma ^ + ^mT^ ^ //ct) fla (^5) 

where t^ is the total probability for the i-site electron with spin a in the m-th orbital 
hopping to the j-site, ttot = ^"i is the total probability for all z-site electrons with 
spin a in all localized orbitals hopping to ?'-site, and tmi = is the probability for the 

(•tot 

i-site electron with spin a in the m-th orbital hopping to the l-th orbital of the j-site, as 
illustrated in Figure [T} and so as in the opposite direction. In this case the numbers of 
outgoing and incoming hopping electrons are identical in each orbital so that the system 
is in an equilibrium state. 

Moreover, when one orbital is not occupied, e.g., = fii^i = 0, the equations for 
the A^-orbital system automatically turn to those for (A^ — 1) orbitals. 

For another lattice, these set of hybridization functions Ama will be calculated by 
an integral over k space for conduction electrons, then be put into the equations to 
calculate the single particle Green's functions in each orbital. 

At the end of this section, we would like to give a discussion of the spin-flip term 
fima'fLfii^'fima and the pair-hoppiug term flma'fimafn^'fu^T which are studied in some 
recent literatures (see, e.g., [331 [31]), but are not studied in this paper. One background 
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(a) Lattice mode! 




(b) Impurity model 



Figure 1. (Color online) Illustration of inter-site multi-orbital hopping in a Hubbard 
type lattice model (a) and hybridizations in a mapped impurity model (b). In 
the Hubbard type lattice model, we only show the hoppings in one direction. If 
neighbouring sites are identical, the hopping in opposite direction will be the same. 

for this is that contemporary LDA-I-DMFT calculations are mostly implemented to 
achieve the correct charge in the strongly localized orbitals and then calculate the 
correct potential for a subsequent LDA calculation so that the spin-flip exchange term 
is less important in this context, which can be recognized from the following. All the 
Green's functions are associated with some kind of physical interactions, e.g., the Green's 
function ^ fifia-fma'i flia ^ labels an inter-orbital fluctuation which corresponds to the 
inter-orbital Coulomb interaction between electrons with identical spin. In turn, the 
spin flip term is associated with a Green's function ^ fjfj/fma'fia', fLa To show the 
influences of the Coulomb interactions and the spin-flip term, we make a comparison: 
for the Green's function ^ fifi^fma] fma ^5 o'^^ has 

{[flflafma, fL] + ) = ^fla, (36) 

while for the spin-flip associated Green's function ^ fj^i fma' fia', fma ^) 

{[fLfma'fla,fLU)=0. (37) 

The above averages (...) are corresponding to their connection to charge which will 
determine the effective Coulomb interaction strengths, the relative position of the two 
Hubbard bands, and the total densities of states in each of the Hubbard bands. One can 
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note that the spin-flip term does not directly connect to charge for a system with the 
above mentioned anti-commutation relations for the operators, whereas the Coulomb 
interactions are associated with charge. It is estimated that the spin-flip associated 
Green's function will mainly contribute to the shape of the Hubbard bands which 
will show more micro structure on top of the Hubbard bands that are determined 
by the charge related Green's functions. Therefore, the Coulomb interactions are 
more important than the spin-flip term so that in this paper we first concentrate 
on the Coulomb interactions, while we have used the mean field approximation to 
treat the inter-orbital Coulomb interactions and have neglected the on-site inter- 
orbital fluctuations. The spin-flip term will be studied in a higher-order decoupling 
scheme where the inter-orbital fluctuations associated with the inter-orbital Coulomb 
interactions will be fully taken into account. As for the pair- hopping term, this requires 
that the electron's starting orbital will be fully occupied and the destination orbital is 
empty. This physical image reflects a rare case and corresponds to the excited states. 
With the same reason as for the spin-flip term, we have neglected to study the pair- 
hopping term in this order of the decoupling scheme. 

3. Results and discussion 

Using the above derived MO-EOM impurity solver, we have investigated a two-orbital 
system with the two-orbital Hubbard model. Here we first present calculated results 
for the case where the inter-orbital hybridization is neglected. In this case, our model 
Hamiltonian in ([T]) will return to the frequently studied Hamiltonian as shown in [21]. 
When the inter-orbital and intra-orbital Coulomb interaction strengths are identical, i.e., 
Uii = U22 = Ui2aa' = Ui2a(j and the Hund's rule coupling constant J = 0, we obtained 
the densities of states (DOS) for the halffilled system, as shown in Figure|2| where /i is the 
chemical potential and we have furthermore assumed that the two occupied orbital levels 
have the relation Ei = E2. We can see that the two orbitals have different broadening 
for the Hubbard bands. With increasing U , the system changes from metallic states to 
insulating states. However, the metal-insulator transitions (MIT) occur at different U , 
not shown in the plot. The critical value of U for the two orbitals are correspondingly 
U = Uci ~ 1.5 and U = Uc2 ~ 2.5 (in units of half the narrow bandwidth). 

In Figure [s] we show the computed results for a non-zero J, J = U/A. We still 
assume Uu = U22, while the inter-orbital Coulomb interaction strength now changes 
to be Ui2acT' = Uii — 2J and Ui2aa = Uu — 3J. Comparing to the identical U case 
shown in Figure [2| one recognizes that with the decrease of the inter-orbital Coulomb 
interaction strength, the effective decreases accordingly, which has influenced the 
MIT in each orbital. In this case, the critical value of U for the two orbitals has 
correspondingly increased to U = Ud ~ 1.9 and U = Uc2 ~ 3.5, for the narrow and 
wide orbitals, respectively. In Figure [2] and Figure [3} we have seen an orbital selective 
Mott transition (OSMT). But due to the neglect of inter-site inter-orbital hopping, here 
the two orbitals are only connected with the inter-orbital Coulomb interaction strengths 
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Figure 2. (Color online) Quasiparticle densities of states for the halfRUed two-orbital 
system on the Bethe lattice in case of neglecting inter-orbital hybridizations. The DOS 
are computed with half bandwidths D2 — 2Di = 2, Hund's rule coupling constant 
J = and at temperature T = 0.01. The red solid line gives the DOS of the narrow 
orbital and the blue dashed line gives that of the wide orbital. 



and occupations. Comparing Figure [2] and Figure |3j when the inter-orbital Coulomb 
interaction strengths decrease, the two Hubbard bands will come closer together so that 
the critical value of U for the two orbitals will increase. Moreover, once the two orbitals 
have identical intra-orbital Coulomb interaction strengths and different bandwidth, the 
OSMT will always occur at the halffilling. 

Next, we consider the two-orbital system taking now non-zero inter-orbital 
hybridizations into account. Doing this, we have calculated the DOS for the 
corresponding two-orbital system at halffilling. First we show the obtained results in 
Figure |4] for the two-orbital system having identical bandwidths for the two orbitals. 
We can observe that the inclusion of the inter-orbital hybridizations has greatly changed 
the DOS. Besides the change of DOS with increasing U, we also give the comparison of 
the DOS at J = and J = U/A. When J is not zero, the effective Coulomb interaction 
strength becomes smaller, which has increased the critical value of U for the Mott 
transitions. 

The quasiparticle weight, which is a characteristic quantity describing the strength 
of the correlation effect, is defined by 

Z = ^=[l-^?55Ml- . (38) 

where is the self-energy, and m and m* the bare and enhanced mass, respectively. 
The quasiparticle weight has been computed as the function of U at halffiling for the two- 
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Figure 3. (Color online) Quasiparticle densities of states for the halffilled two-orbital 
system on the Bethe lattice in case of neglecting inter-orbital hybridizations, computed 
with half bandwidths D2 = 2Di = 2, Hund's rule coupling constant J — U/A, and at 
temperature T = 0.01. The two orbitals are denoted by the dashed and full curves, as 
in Figure [2j 
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Figure 4. (Color online) Quasiparticle densities of states for the two-orbital system 
at the halflilling in case of including the effect of hopping between the orbitals. The 
DOS are computed with half bandwidths Di = D2 — 1 and at T = 0.01. 



orbital system with identical bandwidths. It is shown in Figure |5j The abbreviation 
GA in the legend means that it is calculated with genetic algorithm (see [2B]), and 
"linearmix" means that these data are calculated with the iterative method with linear 
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Figure 5. (Color online) Quasiparticle weight versus U for the two-orbital system 
on the Bethe lattice at halffiUing in the case of including hopping effect between the 
orbitals, computed with half bandwidths Di = D2 = I and at T = 0.01. 

mixing. The value of Z decreases along with increasing the value of U until the critical Uc 
for the Mott transition is reached. However, due to the Lorentzian broadening applied in 
the iterative method with linear mixing, it has given minor different results and smaller 
critical values of U for both J = and J = U/A cases. 

We mention here that all the model calculation results presented in this section, if 
not specified otherwise, are calculated with the genetic algorithm in order to efficiently 
solve the integral equations. Using instead the iterative method with linear mixing 
should have a minor difference on the computed DOS and critical value of U, but would 
nonetheless be less accurate and it will also cost more cpu time. This is so because 
the iterative method with linear mixing has to use a finite Lorentzian broadening 
in calculation of Green's functions on real frequencies. In calculations with the GA 
scheme, this Lorentzian broadening can be set to zero so that the GA scheme will give 
a numerically exact solution of EOMs. For a more detailed explanation, see pO] . 

Next we study the case where we have the halffilled two-orbital system with different 
bandwidths, but identical intra-orbital and inter-orbital Coulomb interaction strengths, 
i.e., J = 0. The obtained quasiparticle DOS is given in Figure [6j The DOS reveals that 
the MIT for the two orbitals occurs nearly simultaneously. Compared with Figure |2| the 
OSMT disappeared. This difference is due to the inclusion of the inter-site inter-orbital 
hopping effect. The critical value of U for the two orbitals is here Uc ~ 2.4. 

When J = U /A the computed DOS are shown in Figure [7j We can see that the 
quasiparticle peak is nearly pinned at the chemical potential for the wide orbital. Near 
U ~ 3.0 the quasiparticle peak for the narrow band has become very small, but it 
still exists. Compared with Figure |3| the OSMT also disappeared under our present 
assumptions and on the Bethe lattice, which can not, however, exclude the OSMT in 
all cases, e.g., if the neighbour sites are different (i.e., the material is a compound) 
or inter-site multi-orbital hoppings are not like what we have assumed in a real system 
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(physical lattices) or for magnetic cases when outer fields exist. But we believe that with 
the introduce of inter-site multi-orbital hopping effect, the OSMT will be suppressed 
in some sense. This phenomena will be investigated in a further study with improved 
techniques beyond the here-made mean field approximation in the impurity solver and 
on other lattices. 

Next we present the quasiparticle weight versus U for the two-orbital system with 
different bandwidths at halffiUing in Figure [8] We can clearly see in Figure [8] that 
the two DOS curves of the wide and narrow orbitals approach zero at the same point, 
Uc = 3.2, which means that the metal-insulator transition for the two orbitals occurs 
practically simultaneously. 

Next we present some comparisons of results computed with the MO-EOM method 
to results achieved with other numerical methods to see how well our MO-EOM impurity 
solver performs. One comparison is visualized in Figure |9} where the DOS obtained with 
our MO-EOM method are compared with those achieved with the NRG method [35] for 
the case that the two orbitals have identical band widths. One can note that both the 
methods agree well with the positions of the Hubbard bands. They do show a difference 
in the critical value of U for the Mott metal-insulator transition, which can however 
be partly explained by the large Lorentzian broadening (which causes a long tail of 
the peak) shown in the NRG data. In our genetic algorithm method this Lorentzian 
broadening can even be set as zero. Considering this broadening effect, our results are 
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Figure 6. (Color online) Quasiparticle densities of states computed for the two- 
orbital system at lialfRUing for the case where the effect of hopping between orbitals 
in included. The results have been calculated with D2 = 2Di = 2 and J = 0, and at 
T = 0.01. The red solid line is for the narrow orbital, and the blue dashed line for the 
wide orbital. 
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Figure 7. (Color online) Quasiparticle densities of states for the two-orbital system 
on the Bethe lattice at halffiUing in case of including the effect of electron hopping 
between the two orbitals. The DOS have been obtained with D2 = 2Di — 2, and 
J = U/A and at T = 0.01. The two orbitals are denoted by the dashed and full curves, 
as in Figure [6] 
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Figure 8. (Color online) Computed quasiparticle weight as a function of U for the 
two-orbital system at halffiUing, obtained with the parameters D2 — 2Di = 2, T = 0.01 
and J = C//4. 



close to those obtained with the NRG method, but may have a shghtly larger critical 
value of U. 

Next we present a comparison of results computed with the MO-EOM method to 



results achieved with the QMC method [31]. This comparison is visualized in Figure 10 
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Figure 9. (Color online) Comparison of quasiparticle densities of states obtained by 
the multi-orbital EOM method with the DOS obtained by the NRG method. The 
system treated is a two-orbital system at halffiUing on the Bethe lattice, with inclusion 
of the effect of electron hopping between orbitals with bandwidth D2 = Di = 1 and 
at temperature T = 0.01, and J — U/A. The DOS shown on the negative ordinate are 
those for the NRG method, those on the positive ordinate for our MO-EOM method. 
The NRG data are taken from [55] . 



though the Hamiltonian we studied here involved the inter-orbital hybridizations, a little 
different to that used in [31]. For visibility purpose, we have plotted the quasiparticle 
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Figure 10. (Color online) Comparison of quasiparticle densities of states obtained by 
the multi-orbital EOM method with DOS obtained by the QMC method. The system 
treated is a two-orbital system at halffiUing on the Bethe lattice, with inclusion of the 
effect of electron hopping between orbitals with bandwidth D2 = 2Di = 2 and at 
temperature T = 1/40, and J = U/A. The DOS shown on the negative ordinate are 
those for the wide orbital, those on the positive ordinate for the narrow orbital. The 
QMC data are taken from [3T] and all parameters are identical to those of |31j . 
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DOS of the wide orbital on the negative y axis and that of the narrow orbital on the 
positive y axis. One can recognize a good agreement between the two DMFT solvers, 
particularly for the lower U values. For the higher U value there is some difference, 
which is related to the fact that the critical value Uc in our calculations is larger than in 
the QMC results. This may however be due to the different Hamiltonians used in the 
simulations. We used a less exact approximation (mean field approximation to inter- 
orbital Coulomb interactions) but a more realistic Hamiltonian, while the QMC method 
used an exact numerical method but a simplified Hamiltonian. Therefore here we can 
not conclude which one is more precise. In our opinion, it is reasonable to obtain a 
larger critical value of U in our MO-EOM method, because we have taken into account 
the hopping between different orbitals of different sites. Moreover, in our treatment we 
have assumed that the probability of electrons hopping from the narrow to the wide 
orbital is ^^^^^^^^ which is much larger than the probability of hopping from one 

narrow to another narrow orbital. If more electrons would be allowed to hop into an 
orbital with identical orbital index, the critical value of U will become smaller. We 
should also note that the mean field approximation made for the inter-orbital Coulomb 
interactions might have narrowed the gap between lower and upper Hubbard band, i.e., 
this would lead to an increase of the critical value of U. Using the present mean field 
approximation to achieve comparable results as the QMC method, the future extension 
beyond the mean field approximation is promising. 

As a next step we investigate the influence of the filling level. To this end we 
have varied the total occupation number in the simulations. We present the computed 
DOS of the two-orbital systems at different occupation numbers in Figure 11, where 
we have taken the parameters that the two orbital levels are identical, i.e., Ei = E2. 
Vertically, the figure shows the changes of the DOS along with the total occupancy. 
When the filling of the orbitals increases, the effective Coulomb interaction is seen to 
become larger. Correspondingly, the DOS can be seen to change from a one or two- 
peak structure to a three peak structure typical of strongly correlated systems, where 
the quasiparticle peak appears at the Fermi level. In this procedure, an orbital selective 
phenomena is observed, i.e., at the beginning the occupation in the narrow orbital is 
less than that in the wide orbital, then changes to equal and then larger than the 
occupation in the wide orbital, finally equal again at the half-filling, ntot = 2. If the 
intra-orbital Coulomb interaction strengths or orbital levels or both are different for 
the two orbitals, the change of occupation and effective Coulomb interaction in the 
two orbitals will be more complicated. In addition we compare in Figure 11 the DOS 
calculated with different Coulomb interaction strengths, U = 2 and U = 3. It can be 
recognized that, when U increases, the two Hubbard bands move farther away from each 
other, as anticipated, and a gap appears in the unoccupied spectrum at, e.g., u — n ^2 
for [/ = 3 and ntot = 1-6. At ntot = 2, the quasiparticle peak is also reducing along 
with the increment of U. When U is large enough, the quasiparticle peak will eventually 
disappear and a gap appears at the chemical potential, as shown in panel (d) of Figure [7| 
and the system will turn to an insulating state. 
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Figure 12 shows the quasiparticle weight versus the total occupation number of the 
two orbitals. Due to the fact that J = U/4 and that the occupation numbers in the two 
orbitals are different, the effective Coulomb interaction for the electrons in the narrow 
orbital is smaller than that for electrons in the wide orbital when ntot is small, this 
causes the narrow orbital to have a larger quasiparticle weight than the wide orbital. 
Yet, when the occupation of the narrow orbital increase, its weight becomes reduced, 
leading to a cross-over point of the quasiparticle weights of the two orbitals. 



4. Summary 

In this paper we have proposed and tested a fast multi-orbital impurity solver for 
the DMFT based on the equations of motion method. In the construction we have 
introduced a mean field approximation for inter-orbital Coulomb interactions appearing 




Figure 11. (Color online) Quasiparticle densities of states computed for the two- 
orbital system on the Bethe lattice for different total fillings. The parameters used are 
D2 = 2Di = 2, T = 0.01, and J = C//4. The red solid line is for the narrow orbital, 
and the blue dashed line for the wide orbital. 
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Figure 12. (Color online) Quasiparticle weights of the narrow and wide orbitals as a 
function of the total filling for the two-orbital system on the Bethe lattice, computed 
with the parameters D2 = 2Di = 2, T = 0.01, and J = U/4. 

with the two and three particle Green's function expression. The proposed impurity 
solver includes inter-orbital hybridizations and it preserves the particle-hole symmetry 
for halffilled states, and it automatically satisfies the sum rule that integrating of the 
density of states over all uj space for each spin and orbital index equals identity. The 
influence of the inter-orbital hopping has been investigated numerically with the MO- 
EOM method. The results obtained with our impurity solver show that the inclusion 
of inter-site inter-orbital hopping does make a difference and indicate that the inter- 
site inter-orbital hopping effect is quite important for multi-orbital systems. It is 
observed that the inclusion of the inter-site inter-orbital hopping effect will strengthen 
the connections between different sites (and hence eliminate the difference between the 
orbitals if the neighbour sites are all identical to this impurity), and has caused the 
disappearance of the OSMT under our present assumptions and on the Bethe lattice. 
This will be tested further in a forthcoming work with treatments beyond the mean field 
approximation and on different lattices. Although our approximations are simplifying in 
some respects, we do believe that, if one wants to study the difference between orbitals 
having different bandwidths, e.g., the OSMT, the inter-site inter-orbital hopping can 
not be neglected. 

Moreover, we note that the decoupling scheme only acts on the on-site inter-orbital 
fluctuation terms of the impurity, and does not influence the spatial fluctuations of 
the lattice model. The spatial fluctuations should be considered in the DMFT self- 
consistency conditions or a future EOM cellular DMFT. 

The precision of the EOM method depends on the decoupling only. A drawback 
of this impurity solver is that the mean field approximation to the on-site inter-orbital 
Coulomb interactions may give rise to the loss of some interesting information when the 
on-site inter-orbital fluctuations are strong. Hence, if one wants to study the inter-orbital 
fluctuations and even phenomena related to higher-order terms (e.g., spin-flip term and 
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pair hopping term), a higher-order decouphng scheme will be needed. Nonetheless, the 
present EOM method can already be applied to most systems. 

The main advantages of the MO-EOM method are, first, that it works directly 
on real frequencies and for arbitrary temperatures, second, it is a fast impurity solver 
for the multi-orbital case. Finally, it can work for an arbitrary number of orbitals. 
Consequently, if one could determine all the parameters between the orbitals properly, 
even all the orbitals of one atom could be treated within MO-EOM method. Therefore 
the developed MO-EOM method is a competitive solver as compared, e.g., to the 
QMC. The multi-orbital EOM solver can be applied to multi-orbital system containing 
several orbital levels, with varying band widths, inter-orbital and intra-orbital Coulomb 
interaction strengths. Also systems with a crystal field having different orbital levels 
can be treated with this MO-EOM method. Besides the model calculation, using this 
impurity solver and the formula that we have given for orbital levels and effective 
Coulomb strength, one can obtain the occupations and then use these to correct the 
charge, orbital levels and remove the double counting, so that it can be implemented in 
DFT-I-DMFT calculations for real correlated materials. 

In the peer review of this paper, our attention was drawn to the paper by P. 
Hansmann et al, who have mentioned the inter-site inter-orbital hoppings in their 
LDA+DMFT work. But these have however not appeared in model calculations to 
discuss the influence of these hoppings on the orbital selective Mott transition. 
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